
rm(list=ls())
library(MASS)
library(mediation)
setwd("C:/Users/dtingley/Dropbox/MyersTingley/analysis/replication")

library(foreign)
library(xtable)
library(mgcv)
library(VGAM)
library(sandwich)
library(Matrix)

load("EmotionsPAReplication.RData")

a<-lm(GD_NegAff~AnxietyVsCont+gender_p+Trust_GSS_pre+ideo7,data=Anxiety )
g<-lm(GD_NegAff~GuiltVsCont+gender_p+Trust_GSS_pre+ideo7,data=Guilt )
f<-lm(GD_NegAff~AngerVsCont+gender_p+Trust_GSS_pre+ideo7,data=Anger )
a1 <- lm(sent ~ AnxietyVsCont+GD_NegAff+gender_p+Trust_GSS_pre+ideo7,  data=Anxiety)
g1 <- lm(sent ~ GuiltVsCont+GD_NegAff+gender_p+Trust_GSS_pre+ideo7,    data=Guilt)
f1 <- lm(sent ~ AngerVsCont+GD_NegAff+gender_p+Trust_GSS_pre+ideo7,    data=Anger)

a2 <- mediate(a,a1,sims=1000,treat="AnxietyVsCont",mediator="GD_NegAff",conf.level=.9)
g2 <- mediate(g,g1,sims=1000,treat="GuiltVsCont",mediator="GD_NegAff",conf.level=.9)
f2 <- mediate(f,f1,sims=1000,treat="AngerVsCont",mediator="GD_NegAff",conf.level=.9)

a<-lm(BasNeg_Fear~AnxietyVsCont+gender_p+Trust_GSS_pre+ideo7,data=Anxiety )
g<-lm(BasNeg_Guilt~GuiltVsCont+gender_p+Trust_GSS_pre+ideo7,data=Guilt )
f<-lm(BasNeg_Hostility~AngerVsCont+gender_p+Trust_GSS_pre+ideo7,data=Anger )
a1 <- lm(sent ~ AnxietyVsCont+BasNeg_Fear+gender_p+Trust_GSS_pre+ideo7,   data=Anxiety)
g1 <- lm(sent ~ GuiltVsCont+BasNeg_Guilt+gender_p+Trust_GSS_pre+ideo7,    data=Guilt)
f1 <- lm(sent ~ AngerVsCont+BasNeg_Hostility+gender_p+Trust_GSS_pre+ideo7,    data=Anger)

a3 <- mediate(a,a1,sims=1000,treat="AnxietyVsCont",mediator="BasNeg_Fear",conf.level=.9)
g3 <- mediate(g,g1,sims=1000,treat="GuiltVsCont",mediator="BasNeg_Guilt",conf.level=.9)
f3 <- mediate(f,f1,sims=1000,treat="AngerVsCont",mediator="BasNeg_Hostility",conf.level=.9)

#save(a2,g2,f2,a3,g3,f3,file="MediationResults.RData")
#load("MediationResults.RData")

pdf("../Combined-NoInt-NoGD.pdf", width=8, height=3, onefile=FALSE, paper="special")
par(mfrow=c(1,3))
xlim<-c(-1.5,1.5)
plot(a3,main="Anxiety Induction/Mediator",xlim=xlim)
plot(f3,main="Anger Induction/Mediator",xlim=xlim )
plot(g3,main="Guilt Induction/Mediator",xlim=xlim)
dev.off()




pdf("../Combined-NoInt-GD.pdf", width=8, height=3, onefile=FALSE, paper="special")
par(mfrow=c(1,3))
#xlim<-range(plot.noint$sent.Anxiety.BasNeg_Fear,plot.noint$Guilt.BasNeg_Guilt)
xlim<-c(-1.5,1.5)
plot(a2,main="Anxiety Induction/GD-Neg",xlim=xlim)
plot(f2,main="Anger Induction/GD-Neg",xlim=xlim)
plot(g2,main="Guilt Induction/GD-Neg",xlim=xlim)
dev.off()


#sensitivity

sens.cont.1<-medsens(a3, rho.by=.1,  sims=1000)
save(sens.cont.1,file="Mediation-Sensitivity-anxiety.RData")

pdf("../Anxiety-Sensitivity-noint.pdf", width=8.5, height=4.5, onefile=FALSE, paper="special")
par(mfrow=c(1,2),mar=c(2, 3, 1.5, 0), mai=c(1,.8,.4,.2), cex.lab=.8, cex.axis=.8 , cex.main=.8, font.main=1)
plot(sens.cont.1, sens.par="rho",main="")
plot(sens.cont.1, sens.par="R2",  r.type="residual", main="",sign.prod=-1)
dev.off()



library(quantreg)
quantile<-seq(.35,.65, by=.05)
med_qr<-matrix(, nrow=length(quantile), ncol=1)
med_up<-matrix(, nrow=length(quantile), ncol=1)
med_low<-matrix(, nrow=length(quantile), ncol=1)
dir_qr<-matrix(, nrow=length(quantile), ncol=1)
dir_up<-matrix(, nrow=length(quantile), ncol=1)
dir_low<-matrix(, nrow=length(quantile), ncol=1)
a<-lm(BasNeg_Fear~AnxietyVsCont+gender_p+Trust_GSS_pre+ideo7,data=Anxiety )
for(i in 1:length(quantile)){
  o <- rq(sent ~ AnxietyVsCont+GD_NegAff +gender_p+Trust_GSS_pre+ideo7, data=Anxiety, tau=quantile[i])
  temp<-mediate(a, o, sims=1000, boot=TRUE, treat="AnxietyVsCont", mediator="GD_NegAff")
  med_qr[i]<-temp$d1
  med_low[i]<-temp$d1.ci[1]
  med_up[i]<-temp$d1.ci[2]
  dir_qr[i]<-temp$z1
  dir_low[i]<-temp$z1.ci[1]
  dir_up[i]<-temp$z1.ci[2]
}
save(med_qr,med_low,med_up,dir_qr,dir_low,dir_up,file="Mediation-Anxiety-QuantileAppendix.RData")

pdf("QuantileRegression-Anxiety-noint.pdf", width=10, height=5, onefile=FALSE, paper="special")
par(mfrow=c(1,2),mai=c(.8,.9,.4,.2) )
g_range<-c(min(med_low),max(med_up))
xlim<-c(min(quantile),max(quantile))
plot(quantile, med_qr, type="n", xlab="Quantile", ylab = "Mediation Effect", xlim=xlim, ylim =g_range)
polygon(x=c(quantile, rev(quantile)), y=c(med_low, rev(med_up)), border=FALSE, col=8, lty=2)
lines(quantile, med_qr, lty=1)
abline(h=0)
mtext("Quantile Mediation Effects", line=1)
g_range<-c(min(dir_low),max(dir_up))
plot(quantile, dir_qr, type="n", xlab="Quantile", ylab = "Direct Effect", xlim=xlim, ylim =g_range)
polygon(x=c(quantile, rev(quantile)), y=c(dir_low, rev(dir_up)), border=FALSE, col=8, lty=2)
lines(quantile, dir_qr, lty=1)
abline(h=0)
mtext("Quantile Direct Effects", line=1)
dev.off()
